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ABSTRACT 

This paper describes a practical approach using a general 
purpose lumped-parameter computer program, GFSSP 
(Generalized Fluid System Simulation Program) for calculating 
flow distribution in a network of micro-channels including 
electro-viscous effects due to the existence of electrical double 
layer (EDL). In this study, an empirical formulation for 
calculating an effective viscosity of ionic solutions based on 
dimensional analysis is described to account for surface charge 
and bulk fluid conductivity, which give rise to electro-viscous 
effect in microfluidics network. Two dimensional slit micro 
flow data was used to determine the model coefficients. 
Geometry effect is then included through a Poiseuille number 
correlation in GFSSP. The bi-power model was used to 
calculate flow distribution of isotropically etched straight 
channel and T-junction microflows involving ionic solutions. 
Performance of the proposed model is assessed against 
experimental test data. 

INTRODUCTION 

Microfluidics channel networks [1,2] are being placed on 
lab-on-a-chip devices to perform multi-step biochemical 
processes, or to perform parallel, arrayed biochemical 
experiments [1, 2], In these microfluidic channel networks, 
typical length-to-width aspect ratios are very high (usually in 
thousands). Traditional CFD (computational fluid dynamics) 
methods are not practical for microfluidics-network design. 
System level (or lumped parameter models) tools [3, 4, 5] were 
recently developed for modeling electro-osmotic pumping, 
diffusion, lab-on-a-chip layout and design. When considering 
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liquid (particularly electrolytic) flows in micron-sized channels, 
interfacial effects due to the existence of the electrical double 
layer (EDL) cannot be ignored. Surface charge occurs when 
solid surfaces are in contact with aqueous solutions due to the 
asymmetric dipoles of the molecules that reside at the interface. 
In microfluidic networks, the EDL induced electrokinetic 
effects, namely, electro-osmosis and streaming potential are 
well known [6], In the absence of external applied electrical 
potential, an important electrokinetic effect for pure pressure 
driven micro flows is the so-called electro-viscous retardation 
effect. When a pressure gradient is applied to fluid flow 
through a micro channel, a streaming electric field is generated 
due to the EDL. This streaming potential is counter directional 
to the fluid and causes flow retardation. This phenomenon 
manifests itself in an apparent viscosity, which is larger than 
the actual fluid viscosity according to several theoretical and 
experimental studies [7-11]. The purpose of this study is to 
develop a simplified electrokinetic model based on dimensional 
analysis to account for the electroviscous effect in pressure 
driven microfluidic network. The model is validated against 
experimental studies. 

MODELING APPROACH 

To account for the electro-viscous effect, a body force 
term, resulting from the existence of surface charge and the 
resulting EDL, is added to the fluid momentum equations 
(neglecting gravity force in typical microfluidic networks), 
written as: 
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The last term of the equation (1) represents the electrokinetic 
effect. For streaming potential induced by the applied pressure. 


K E , = -U,p e , ( 2 ) 

and the local electric charge density, p e , is governed by the 
Poisson equation 

V> = -Bs - . (3) 

££ 0 

In the above equations, p is the fluid density, u ; the Cartesian 
velocity components, p the pressure, x tJ the stress tensor, s the 

fluid permittivity, and A 0 the bulk fluid conductivity. The 
charge density is expressed as p e = e(Z + n + Z n_) and the 


ionic number concentration, n ± , can be obtained by solving the 


Nemst-Planck equation. Under equilibrium situation, the ion 
distribution can be described by the exponential function; 
equation (3) becomes the well-known Poisson-Boltzmann 
equation: 
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where L is the length of the channel, A P is the pressure 
difference, £ 0 is the zeta potential, and V is the fluid kinematic 
viscosity. A bi-power law form was used such that 


^=i+a(n,) b (n ! )° (6) 
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To determine the model coefficient, the slit microflow 
experiments of Ren et al. [11] were used to determine the 
constants a, b, and c in the functional relationship for the 
effective viscosity, yielding the correlation: 


-^ = l + 0.856(n i ) 1 ’ 1 °(n 2 ) 1 ' 72 ' (7) 

The slit flow of [11] was assumed to be two dimensional. To 
account for various channel cross section geometries 
encountered in microfluidic networks, the approach utilized in 
GFSSP (Generalized Fluid System Simulation Program, [14]) 
was used through the Poiseuille number correlation. This 
procedure consists of the following steps: 

1 . Estimate the hydraulic diameter of the cross-section: 

D h = (4)(Area) / Perimeter 


The characteristic EDL thickness can be defined based on the 
inverse of the Debye-Huckel parameter K , and 
k 2 =(2Z 2 e 2 n 0 )/(££ 0 k B T). Assuming Newtonian fluids, the 

net effect of the EDL (last term of equation (1)) would produce 
extra resistance force to the fluid flow. Previous theoretical 
studies strive to derive analytical form of an apparent or 
augmented viscosity. Most assumptions involved fully 
developed flow field, constant fluid properties, and most 
importantly, the Poisson-Boltzmann equation (either the linear 
or non-linear form) for EDL. With the Boltzmann assumption, 
in which the ionic concentration follows an exponential 
distribution, the analytical forms of apparent viscosity usually 
have explicit Debye-Huckel parameter ( K ) dependency (see 
[7-10]). Since K is proportional to the square root of ionic 
molar concentration, formulations based on the Boltzmann 
assumption would lead to singularity for very diluted 
electrolyte solutions (see the discussion in [10] and chapter 3 of 
[7]). To avoid this anomaly, the Nemst equation was used to 
obtain the ionic concentration distribution [12] and a close- 
form formulation cannot be obtained. 

In this study, a different approach was taken. We proposed 
an empirical form of the apparent viscosity formulation using 
surface charge and fluid’s bulk conductivity as two extra 
dependent variables to account for the electro-viscous effect 
based on dimensional analysis [13]. The two dimensionless 
groups based on the above-mentioned parameters are: 


2. Estimate the Poiseuille number for a particular cross-section. 
The Poiseuille number can be expressed as a polynomial 
function of aspect ratio. Coefficients for different geometries 
can be found in [14]. 
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3. Calculate the friction factor for a non-circular pipe: 

_ 4 Po * f or L am j nar f] ow (Re<2300) (9) 

Re 

GFSSP employs a finite volume formulation of mass and 
momentum conservation equations in a network consisting of 
nodes and branches. Mass conservation equation is solved for 
pressures at the nodes while the momentum conservation 
equation is solved at the branches to calculate flowrate. For 
given pressure drop across various ports of the microfluidic 
network, the GFSSP program calculates fluid velocities in each 
branch and pressure at each node. Further detailed numerical 
method utilized in the GFSSP program can be found in 
Majumdar et al. [14], 

EXPERIMENTAL PROCEDURE 
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Microfluidics experimental study was carried out on a 
quartz chip using the Application Development Unit (ADU). 
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The quartz chip consists of two joined layers. The bottom layer 
contains seven isotropically etched channels with an identical 
height and width of 15 and 50 pm, respectively. The top layer 
has 25 wells that align with the endpoints of the microchannels. 
The resulting chip thickness is 1 .462 mm and the isotropically 
etched channel cross-section can best be described as bowl- 
shaped. Figure 1 shows the various microchannel 
configurations. In this study, only the straight channel (label 
number 4) experiment was performed. Housed by the ADU, the 
chip is placed on a frame-like platform which resides above 
three, face-up, movable microscope lenses. A head gasket 
directly above the platform provides a seal when lowered onto 
the chip. Specific pressures can then be applied to any of the 
25 wells. The ADU has a pressure range of -5 to 5 psi. The 
setup is shown in Figure 2. 

A data acquisition and GUI control system (LabView) is 
used to manage the ADU. Well pressure, lens magnification, 
focus and movement, and light intensity can be adjusted by the 
control system. The magnified microchannel image is 
displayed on an adjacent monitor, and fluid movement can be 
observed during operation. Hash marks of known length along 
certain channels were used to indicate start and stop times for 
subsequent velocity calculations. 

Two different concentrations of aqueous KC1 solutions 
were used as testing liquids: HE 2 and 10" 4 M KCL. The 
addition of dye was necessary for visualization of fluid 
movement, particularly the interface, along the channel. 
Aliquots of testing liquid were added using a 20 pL pipette to 
the wells at either end of the channel to minimize capillary 
forces. The chip was then placed in the ADU and specific 
pressures were applied to one well to initiate fluid movement to 
the lower pressure well. To create the interface, fluid at the 
higher pressure well contained the dye while the fluid at the 
lower pressure did not. A stop watch was used to measure the 
interface travel time for a predetermined distance. Additional 
travel distance was allowed for the fluid to reach steady state. 
After passing the end marker, the direction of fluid flow was 
reversed by lowering the pressure on the well. 

Pressure readings from the ADU were used to determine 
the pressure differential. Multiple runs at the same pressure 
were carried out on the same channel. The channel was 
cleaned with 1 M NaOH and then DI water before changing the 
testing liquid. The process was repeated for a range of 
pressures for each of the testing liquids. Only the straight 
channel results are presented in this paper. 


RESULTS & DISCUSSION 

The proposed model (Eq.7) was implemented into the 
system level computer program GFSSP for predicting the 
straight microchannel flow (# 4 in figure 1). Table 1 shows the 
predicted results versus the experimental data. As seen from the 
table, the reduction of flow velocity due to electro-viscous 
effect can be observed experimentally. The degree of flow 
retardation is quite small for the pressure drop range used in 
this experiment. The model predictions also show flow 
reduction as expected. However, the degree of flow retardation 
was overpredicted. The data indicated about 10 percent 
reduction in flow velocity (as indicated by lower Reynolds 
number) from 10' 2 M to 10~ 4 M KC1 solution flows. The model 


results show larger, up to 20 percent, reduction from no-EDL- 
effect to EDL-effect (10‘ 4 M KC1) flows. In this study, the 
flowfield generated by 10' 2 M KC1 solution was assumed to be 
unaffected by electro-viscous effect following the suggestion of 
Ren et al [11]. Due to the difference in channel shape, the 10" 2 
M KC1 flow may exhibit electro-viscous effect in this study. 
This, however, was not modeled here. For the 10' 2 M KC1 case, 
laminar simulation results (without EDL effect) were used for 
comparison. The pressure drops used in this study were also 
much lower than those used in the slit-flow experiment of Ren 
et al. [11]. At higher pressures drops, the current flow 
visualization technique used for flow velocity measurement 
becomes unreliable. The model did show consistent trend in 
this straight channel case. 

Next, we turned to the application of the current model to 
the T-junction micro flow on the Caliper N145 chip [15], 
Figure 3 shows the layout and dimensions of micro-fluidic chip 
considered in the present study. The glass chip was micro 
fabricated by isotropic etching with about 60 micrometer in 
width and 12.4 micrometer in depth. The corresponding GFSSP 
model is shown in Figure 4. By adjusting the pressure 
difference among well number 1, 8 and 4, different flow fields 
across the T-junction can be set up. Dyes were introduced in 
well number 1 and well number 8 for flow velocity 
measurement [5] using vitalization technique. To accommodate 
flow visualization, pressure levels at well number 1, 8 and 4 
were always different with well number 4 being the minimum. 
For micro-channel with this dimension, the electro-viscous 
effect would be significant for water and low concentration 
ionic solution fluid flows. In addition to pure water, aqueous 
KC1 solutions of concentration 10' 4 and 10' 2 M were used as the 
testing liquids. Following the previous studies [7, 11], the 
10‘ 2 M experiment was taken to be the case without EDL 
effects when compared to the numerical results. Figure 5 shows 
the comparison between numerical prediction and test data. 
The velocity in branch 34 is plotted with the maximum pressure 
differential between inlet and outlet. Experimental data has 
been compared with GFSSP prediction. The variation between 
the 10" 2 and 10' 4 MKC1 velocities increases as the pressure 
difference increases. At higher AP values, conventional laminar 
result overestimates the velocity of the 10" 4 M KC1 solution, 
which is considered to be affected by the EDL. When the 
electro-viscous effect is taken into consideration, the results 
show good comparison against experimental data for 10' 4 M 
KC1 with no more than a 5% difference. The electroviscous 
effect accounts for a decrease in velocity of up to 17% for the 
experimental data. Therefore, while the 10" 2 M KC1 solution 
was not affected by the EDL, the 10' 4 M KC1 solution required 
GFSSP to consider the electroviscous effect for better 
predictions. 


CONCLUSIONS 

An empirical model suitable for system-level lumped 
parameter analysis was proposed to account for electro-viscous 
effects in microfluidic network. The electro-viscous effects 
arise from the electrical double layer (EDL) naturally formed 
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between flowing ionic solutions and solid walls. Based on 
dimensional analysis, the current model includes surface charge 
and fluid bulk conductivity as two extra parameters, in addition 
to conventional parameters used in pressure-driven laminar 
flows. The proposed bi-power form was used to calculate an 
augmented viscosity and a Poiseuille number correlation was 
used for channel geometry effect. The model was tested against 
experimental data of flow distribution in isotropically etched 
straight channel flow and T-junction microflow involving ionic 
solutions. For the straight channel microflow, the flow rate 
retardation of 10% was observed experimentally. The model, 
however, overpredicts the retardation in straight channel flow. 
Explanation of this discrepancy was discussed. In T-junction 
microflow, as much as 17% flow retardation was observed in 
high pressure drop case, and the model predictions compare 
favorably for this case. Further improvement on flow velocity 
measurement is required. Testing of the model in more 
complicated networks and extension of the model to include 
entrance effect are the subjects of on-going studies. 
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Figure 2. ADU platform 



Figure 3. Micro-fluidic chip considered for model verification 
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Figure 4. GFSSP model of flow network in chip with micro-channels 
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* KCI 10-2M exp. 
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x pred.-No EDL 
a pred.-with EDL 


Figure 5. Comparison between current model predictions and test data 



